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Abstract 


The radiation field over a broken stratocumulus cloud deck is simulated by the 
Monte Carlo method. We conducted four experiments to investigate the main factor for 
the observed shortwave reflectivity over the FIRE flight 2 leg 5, in which reflectivity 
decreases almost linearly from the cloud center to cloud edge while the cloud top height 
and the brightness temperature remain almost constant through out the clouds. From our 
results, the geometry effect, however, did not contribute significantly to what has been 
observed. We found that the variation of the volume extinction coefficient as a function 
of its relative position in the cloud affects the reflectivity efficiently. Additional check of 
the brightness temperature of each experiment also confirms this conclusion. 

The cloud microphysical data showed some interesting features. We found that 
the cloud droplet spectrum is nearly log-normal distributed when the clouds were solid. 
However, whether the shift of cloud droplet spectrum toward the larger end is not 
certain. The decrease of number density from cloud center to cloud edges seems to have 
more significant effects on the optical properties. 


Chapter 1 


INTRODUCTION 


The recent concern about the possible effects of anthropogenic pollutants, such as 
C0 2 , CH 4 , CFC, etc., on the climate change has urged the study of feedback mechanisms 
within our climate system. Cloud feedback mechanisms, which include changes in cloud 
type, coverage, and cloud droplet size induced by the change of some other climatic 
factors, are among the least understood feedback mechanisms. Yet, the importance of 
the effects of clouds on the radiation budget has been acknowledged. The magnitude of 
the global shortwave, longwave and net radiative forcing are all one order greater than 
that of C0 2 doubling (Ramanathan et al. ,1989). This implies that small changes in the 
cloud radiative forcing could result in a significant alternation of the Earth's radiation 
budget. 

From a radiation budget point of view, stratus and stratocumulus clouds have 
long been recognized for their importance to the climate system because of their large- 
area coverage and longevity. Both of the stratus and the stratocumulus clouds are 
typically found on the eastern rim of sub-tropical high pressure systems: for instance, the 
west coast of California in the summer, southern Peru, northern Chile, and south west 
Africa. Yet, the spatial scales of stratocumulus and stratus clouds are different: an 
individual stratocumulus cloud ranges around 10 km while the stratus deck could 
extend to 10 6 km 2 in area. 

However, one of the defects in all the present general circulation models (GCM) 
is their inability to calculate radiative transfer through the broken cloud field, which is a 
feature of stratocumulus clouds. Most of the GCMs assume that a homogeneous cloud 


deck covers the entire grid point once the grid point is saturated. It is postulated that the 
radiation budget over a broken cloud field would differ significantly from its 
homogenous counterpart. 

Many efforts have been devoted to the investigation of radiative transfer through 
inhomogeneous clouds. Blerkom (1971) studied how the striate cloud-top structure 
which could result from convective motion in the stratiform clouds affects the shortwave 
radiation field. He found that the reflected radiance depends on both the degree of 
striation and the incident angle. 

McKee and Cox (1974) investigated the reflected irradiance of a cubic-shaped 
cloud and compared it with that of horizontally semi-infinite clouds. Their analysis 
showed that the incident solar flux on the vertical sides of cubic clouds has a significant 
effect on the reflected irradiance. 

Welch and Wielicki (1984) surveyed the effects of cloud shapes and cloud 
alignment patterns on radiation. They studied the bulk reflectivity at visible wavelengths 
as a function of the incident solar zenith angle and fractional area coverage for four 
different shapes of clouds aligned in a linear or hexagonal array. Their conclusion was 
that plane-parallel calculations are not satisfactory at most values of cloud cover. Hence, 
the importance of cloud shapes and cloud alignment patterns to the radiation field is 
assured. 

Yet, despite the inclusion of cloud geometry effects, these previous studies 
assume the optical properties, especially the visible volume extinction coefficient, remain 
constant in clouds. However, it is conjectured that turbulent motions in stratocumulus 
clouds, rising motion in the cloud center and sinking motion at the cloud edge, could 
result in horizontal and vertical inhomogeneity of cloud microphysics properties. Thus, 
the volume extinction coefficient might vary horizontally and vertically in stratocumulus 


clouds. 


Furthermore, these previous studies mainly focus on the study of bulk reflectivity. 
This is understandable because, from a climate point of view, it is the bulk reflectivity of 
clouds that relates directly to the radiation budget. Besides, no in situ measurements of 
cloud point reflectivity, which is defined as the upward flux at a point above the cloud 
dividing by its associated downward flux, were available then. 

The Marine Stratocumulus Intensive Field Observations (IFO) phase of the First 
International Satellite Cloud Climatology Project Regional Experiment (FIRE) in the 
summer of 1987 collected precise and interesting data of stratocumulus clouds. Among 
these data, we are primarily interested in the change of the cloud point reflectivity along 
the flight path detected by the NCAR Electra on June 30. The cloud point reflectivity 
changes drastically from cloud edges to cloud centers in a wave-like form with troughs at 
cloud edges and ridges at cloud centers. In contrast, the cloud-top height remains almost 
constant except a shaip decrease of height at cloud edges. 

We can not help but ask the question: is this a result of the shape of the cloud? 
Or are some other factors associated with this phenomenon? Note, that here the bulk 
reflectivity is not the concern, but rather the cloud point reflectivity as a function in 
space. 

The Monte Carlo method is the only known computational method capable of 
simulating radiative transfer through media in which complicated inhomogeneity is 
present. The reason is that, unlike the other computational methods for atmospheric 
optics, the Monte Carlo method simulates the radiation field by tracking photons, which 
gives the Monte Carlo method flexibility in dealing with complicated inhomogeneity in 
the media. Thus, for cases like simulating the cloud point reflectivity over broken cloud 
field with the change of cloud optical properties considered, the Monte Carlo method is 
indeed the only choice. As a matter of fact, it terms out that tracking photons is the easy 
part of the simulation whereas specifying the inhomogeneity of the cloud field and 
discussing the statistics of results are more challenging. 


Chapter 2 briefly describes the basic concepts and the equation of radiative 
transfer in atmospheric science. The theory, procedures, and statistics of the Monte 
Carlo method are also illustrated for their application to simulations of the radiation field 
of both homogeneous and inhomogeneous clouds. 

Chapter 3 introduces the instruments used in the IFO phase of FIRE, and also the 
general cloud conditions when the data was collected. Some basic arithmetical 
calculations implied in the data are shown. The data of the cloud microphysics 
properties and the cloud radiative filed arc discussed. 

Chapter 4, first of all, briefly describes important concepts about the Mie 
scattering code, the hemispheric mean technique, and the flux extrapolation scheme. 
These three techniques help the fulfillment of our simulations of the cloud radiation field 
and the comparisons between simulation results and observational data. Second, it lists 
all the constants and variables in the simulations. Detailed display and discussions of our 
simulation results are then demonstrated. A concise conclusion of this study is presented 
in Chapter 5. 


Chapter 2 


Computation Method 


In the field of atmospheric radiation, the equation of transfer, which can be solved 
analytically only in some special cases, is noted for its mathematical complexity. Thus, 
efforts have been made on the development of computational procedures. 

Among these techniques, the Monte Carlo method is most capable of simulating 
the radiative field of horizontally inhomogeneous clouds, especially when the cloud 
geometry and the horizontal variation of the volume extinction coefficient are both 
considered. Therefore, the Monte Carlo method is employed to simulate the radiation 
field of stratocumulus clouds in this study. Chapter 2 will briefly review the basic 
concepts of the equation of transfer, and the theory, procedure, and statistics of the 
Monte Carlo method. 


2.1 The Equation of Transfer 

Intensity (radiance), in units of energy per area per time per frequency and per 
steradian, is a basic radiometric quantity in the radiative transfer process. Even though 
intensity is a scalar, a directionality is implied in its definition. Figure (2.1) is a 
schematic diagram of the vectors related to intensity. I x (r,co) denotes the intensity of a 

pencil of light at position r with a wavelength X pointing with a direction co (co is an unit 

vector) at point A. However, as the pencil of light passes through the medium from A to 
B, its intensity changes to I x (r + dr,(o), where dr = cods and ds is the infinitesimal 

distance between point A and point B.. 


z 



Figure (2.1) A schematic diagram for a beam of light,, intensity I^(r,ca) at point 
A, passes through a medium from point A to point B and its associated directionality in an 
arbitrary 3-dimensional coordinate. 




Compare the intensity at point A with the intensity at point B. One will find out 
that the difference between I x (r+dr,(o) and I^(r,(o), 

dl^(r,©) = I^(r + dr,G))-Ij l (r,co) (2.1), 

is a result of loss and gain of intensity from three different mechanisms. 

As the light beam originally with a value I^fr,©) traverses through the medium 

from A to B, it is attenuated because a portion of the light beam is absorbed or scattered 
out into another direction by the medium. The monochromatic volume extinction 
coefficient b ext ^(r), which is the sum of the monochromatic volume scattering 

coefficient b scat ^(r) and the monochromatic volume absorption coefficient b abs ^(r), 
characterizes the amount of attenuation. Note b exa (r), b scat X (r), and b abs>x (r) are all 

in units of per length. The single scattering albedo GJ 0 is defined as ■ scat *^ — . Hence, if 

W(r) 

only the loss term is considered, 

dl^r.GJ) = -b extX (r)I x (r,co)ds (2.2). 

The gain mechanisms add complexity to this problem. We will introduce the 
concept of a light beam as a packet of photons in order to simplify our discussion. 

In 1905, Einstein proposed the photon theory: basically, energy of 
electromagnetic wave behaves as if it is concentrated in photons instead of spreading out 
uniformly over the wave fronts. Based on this theory, a beam of light can be seen as a 
packet of photons. Apply this notion to radiative transfer, the interaction between the 
electromagnetic wave with media which the wave passes through becomes photons 
colliding with particles of the media. 

Previously, ds is denoted mathematically as a scalar with an infinitesimal value. 
From a physical point of view, we can define ds as a distance much shorter than that of 


the mean free path of collisions between photons and particles of the medium. Thus, 
within ds, the probability that more than one collision takes place is very small. Then we 
can simplify our discussion about multiple scattering which occurs within the path ds. 

First, photons which before the collision with particles of the medium propagate 
in the direction co' could be scattered into to after the collision. P(r,(o\co), the phase 
function, is the probability of this event. Note the phase function is also a function of 
position. Therefore, within ds, this gain of intensity in direction co can be expressed as 
bscaJ P(r,co',co)dco'. 

Second, emission of light by the particles in the medium could also contribute to 
the gain of intensity. This term could be written as <J>(A,,r,G)). 

Combine the three mechanisms above, we can write the equation of transfer as: 

2 it 

co-grad I^(r,co) = -b ext> x(r)Ix(r,CD)ds + bscat,X.( r ) Jl A (*\© ) PC 1 "’ 01 ,co)dco 


+ G>(X,r,co) (2.3) 

For homogeneous media, in which b ext j, bscat,l> P and d> are not functions of r. 
Equation (2.3) may be solved analytically. However, solving this integral-differential 
equation analytically thus to obtain the radiation field is almost impossible for 
inhomogeneous media. Therefore, many efforts have been devoted to solving the 
equation computationally. However, most of the computational procedures are limited to 
the plane-parallel atmosphere approximation, in which all the variables in equation (2.3) 
only varies with height. 

Horizontal homogeneity might be a good assumption for the simulation of the 
radiation field of stratus, which often covers area of 10 6 km 2 . However, this assumption 
is not suitable for stratocumulus clouds, in which cell structures are often observed. The 


next section elaborates on the theory, procedure and statistics of Monte Carlo method and 
explains how it simulates the radiation field of stratocumulus clouds. 


22 . The Monte Carlo Method for An Infinite. H omogeneous. Non- Absorbing Cloud 
Deck 


Instead of solving the equation of transfer directly, the Monte Carlo method 
approaches the problem of radiative transfer by modeling stochastic collisions between 
photons and particles. By direct simulation of photon trajectories in the cloud, the 
reflectivity or transmissivity of the cloud becomes one of determining what percent of 
the incident photons are reflected or transmitted under the premise that a sufficient 
number of incident photons has been input. 

For simulating an infinite, homogeneous, non-absorbing cloud deck, in which the 
monochromatic volume extinction coefficient (abbreviated as the extinction coefficient in 
the following text), phase function are constant through out the whole cloud deck, the 
Monte Carlo method could be less efficient than other computational method because it 
costs computer time for optically thick clouds. Nevertheless, a close look at the way how 
the Monte Carlo method simulates light going through an infinite homogeneous cloud 
deck provides insight to the method. The theory, computational algorithm and statistics 
are discussed in the following sub-sections. 

2 21 Theory 

The Monte Carlo method takes full advantages of the photon theory in solving the 
problem of radiative transfer. Since a beam of light behaves like a packet of photons, we 
can simulate photon collision processes, namely, photons colliding with cloud droplets in 
this case, and obtain the radiation field. The whole simulation is only a photon tracing 


scheme in which only simple geometric computations are necessary. The idea of the 
Monte Carlo method is straight forward; however, in order to simulate collision 
processes faithfully, specifying the optical properties of cloud droplets becomes decisive 
to simulation results. The optical properties necessary to be specified for a non- 
absorbing cloud deck, in which the single scattering albedo is equal to 1, are the mean 
free path, the probability density function for the free path, the asymmetry factor, and the 
probability density function for the scattering angle. 

l 0 , the mean free path of collisions between a photon and a cloud droplet, is 
defined as 1 / b exU . b cxU stands for the extinction coefficient, which is the integral of 

the volume extinction cross section of all cloud droplets over a unit volume. In another 
word, when paving all the cloud droplets of a unit volume onto a unit surface, b exl * v 

equals the ratio of integral of cross section to the unit area. Thus, the inverse of b ext 

stands for the height of a unit-area-based volume for which the integral of cross section is 
equal to one unit area. Hence, as a mean, when a photon travels downward from the top 
of this volume, it will collide with one cloud droplet before the photon leaves the 
volume. That is why l 0 is defined as such. 

The Poisson process suits to describe collisions between a single photon and 
cloud droplets, because these collisions possess the necessary and sufficient properties of 
the Poisson process. (1.) They don’t have memory to previous collision events. (2.) The 
probability that a single collision would occur in a small region is proportional to the size 
of the region and does not depend on the number of collisions outside this region. (3.) 
The probability that more than one collision would take place in such a small region is 
negligible. 

The probability distribution of the Poisson random variable x, representing the 
number of outcomes occurring in a given interval or a specified region denoted by t, is 
given by 



p(x,vt) = 


(vt)*exp(-vt) 


x = 0,1,2,. 


(2.4), 


X! 


where v may be interpreted as the mean number of events per unit time or per unit 
volume. Applying this formula to collision processes, we will get equation (2.5) 


p(x ±) = J° — SL. x = 0,1,2,... (2.5), 

l 0 x! 

in which x is the number of collisions within the region, t is the size of an a region. 

The free path / is defines as the distance between two successive collisions. For 
an arbitrary distance t, the probability that the free path / exceeds t, p c (l > t), equals the 

probability that no collisions occur within t, p c (l >t) = p(0,y-) = exp(-y-). Therefore, 

l o l o 

the cumulative distribution function p c (0 < / < t) becomes 

p c (0<,l<,t) = 1- p c (/>t) = 1 - exp(-~) (2.6). 

l o 

We may differentiate equation (2.6) with respect to t to obtain the probability density 
function for free path: 


p(/) = ~‘Cxp(-*— ) (2.7) 

‘0 l 0 

The asymmetry factor is defined as g = cos0 =— — J 0 2 *1 0 X P(cos0)cos©sin0 d© d(j) 

4jc 

, where P(cos0) is the phase function and 0 is the zenith angle and <|> is the azimuth 
angle of the scattering vector (0,<J>). The scattering vector is defined relative to the 
direction in which the photon propagates before a collision. 0 is between 0 and nf 2 for 
forward scattering, between k/2 and 7t for backward scattering. <f> is 0 when the scattering 



vector points at the positive x direction of our arbitrarily defined (x,y,z) coordinate, and 
<j) increases counterclockwise to 271. 

The Henyey-Greenstein phase function is used in this research: 

P (cos©) = (1 - g 2 ) / (1 + g 2 - 2g cos0) 3/2 (2.8) 


2.2.2 Comp utational Procedure 

Figure (2.2) is the flow chart of the procedure of the Monte Carlo method for an 
infinite, homogeneous, and non-absorbing cloud deck. From the chart, the essence of the 
Monte Carlo method is no more than tracking incident photons. The only condition is 
that we must know g, l 0 , (Do first in order to proceed tracking. 

(D 0 determines the probability that a photon is absorbed or scattered when it 
collides with a cloud droplet. For visible radiation through clean water clouds, U5 0 ~ 
0.9999. Thus, for optically thin water clouds, m 0 = 1 is a good approximation. 

From the previous section, we know that the probability density function for the 
free path, /, is an exponential distribution (see equation (2.7)). Since the integral of 
equation (2.7) from / = 0 to / -**> is equal to 1, we can define a one-to-one relation for a 
random number Ry, whose value is uniformly distributed between zero and unity, and a 
free path as follows: 

/ = / 0 In (1 / Ry) (2.9). 


In other words. 




no 


the photon is 
transmitted 


bottom 



top 


the photon is reflected 


Figure (2.2) The flow chart for simulating the reflectivity and transmissivity for a 
infinite, homogeneous, non-absorbing cloud deck, g, l 0 , G3 0 , and b cxU ^ are known. 















Free paths are thus generated by the use of equation (2.9). 

The phase function, which is the probability density function for the scattering 

process (— j 0 2 “/ 0 * P(cos0)sin0 d0 d<f> = 1), is then divided into 1000 equal probability 
4tc 

J 

bins. Each bin has the same area in the probability density function versus solid angle 
diagram. For each bin, we can define a mean angle 0 m : 

— ^OS { — [COS0 (at t}, c left margin of the bin) COs0^ at jj, c right margin of the bin)! ) (2* ^)- 

Thus, when a bin is randomly picked up, 0 m corresponded to the bin is set to be the 0 
angle of the scattering vector. 0 is much easier to obtain: <f> = 2 7 t R^, where is a 

random number between zero and unity. 

By tracking photons, we then can count the number of photons that exit from the 
top of cloud, N out]top . The cloud bulk reflectivity, Ref b , whose definition is the fraction 
of incident radiation reflected, is simply calculated as 

Refb = N out , top /N^ (2.12), 

where is the number of incident photons. 

One important thing worthy of mentioning is equation (2.12) does not contradict 
to the common definition of reflectivity 

Ref b = Ff dd 

ttop /F^id »top (2.13), 

where F T cldtop is the upward flux density (irradiance) at cloud top, and F^ cld top is the 
downward flux density at the cloud top. 



2.2.3 Statistics 


Another question arises naturally after equation (2.12). How many incident 
photons do we need to input in order to get a satisfying reflectivity? Or how many 
incident photons are necessary to make the standard deviation of reflectivity small 
enough? To answer this question, we have to move our focus from the random walk of 
photons inside the cloud to the definition of reflectivity. 

In the case of infinite homogeneous cloud, cloud reflectivity at any point, which 
is constant all over the cloud, has the same value as the cloud bulk reflectivity. 
Regardless of the complicated random walk trajectories of photons, reflectivity is the 
probability that an incident photon is reflected. From this point of view, an incident 
photon being either reflected or transmitted is similar to toss a coin whose chance for 
head or tail is uneven; say, head is for a photon being reflected. Thus the Gaussian 
approximation to the binomial distribution is applicable to describe its statistics (Welch 
and Wielicki, 1984; Lenoble, 1985): 

a = [Ref b (1 - Ref^ / N^] 1/2 (2.14) 

, where <5 is the standard deviation. The largest value of o for the same number of 
incident photons occurs at R = 0.5. When R = 0.5 and <5 = 0.01, is about 2500. 

Table (2.1) lists the reflectivity and transmissivity computed by our Monte Carlo 
code in comparison with those computed by the doubling method. Only 2500 photons 
were input for each Monte Carlo calculation. Our results show good agreement with the 
doubling method especially when the optical depth is more than 0.25. Note that we have 
listed both the value and the percentage of discrepancies. Special caution is necessary for 
those cases with small reflectivity; the percentage of differences may be drastically large, 
even though the absolute value for differences are not off much. However, the fact that 



Table (2.1) A comparison of reflectivity calculated by the Monte Carlo method 
and the doubling method. The difference is calculated by reflectivity (Monte Carlo 
method) minus reflectivity (Doubling method). The difference in percentage form is 
calculated by dividing the difference of the two methods over the reflectivity calculated by 
the doubling method. |i 0 is equal to the cosine of the incident zenith angle. In all of these 
calculations, G3 0 = 1 and g = 0.75. 

























































the larger discrepancies occur at all the small optical depths (when the optical depth is 
0.25) or reflectivity is less than 0.1 instead of when the reflectivity is around 0.5 (for 
example, when the optical depth is 4 and the cosine of the incident zenith angle is 0.5, or 
when the optical depth is 1 and the cosine if the incident zenith angle is 0 . 1 ) suggests 
some deficiency in the application of the Gaussian statistics to the Monte Carlo method. 

.2.3 The Monte Carlo Method for Horizontally Inho mogeneous. Non-absorbing 
Clouds 

2.3.1 Comp utational Procedure 

Our use of the term "horizontally inhomogeneous cloud" primarily denotes that 
b e xt,x varies horizontally or / and cloud geometry needs to be considered. The logic of 

the procedure for horizontally inhomogeneous clouds is pretty much the same as that for 

the homogeneous cloud deck. The differences are described below. 

When b ext ^ is not constant, the mean free path is not constant either. For 

simplicity's sake, consider first a photon traveling from medium A to medium B (Figure 
2.3. a), where b ext is the extinction coefficient for medium A, and b ext for medium 

B. Both b cxt x,a and b ext ^3 are constant Following equation (2.8), a somewhat more 
complicated expression for the photon free path / can be determined: 

/ = OF = OM + MF 

= OM (1 - bexa ' ~) + — - — In (— ), (2.15) 

bext,^3 b cx t,A,B 

where Rj is a random number between zero and unity. 


preceding page blank not filmed 




Figure (2.3. a) A schematic diagram of a free path: O is the position of the (i-l) 
collision; M is the interface point of medium A and medium B in the free path; F is the 
position of the i collision. F is the i collision location if the free path is not calculated 
discretely. 



Figure (2.3. b) A schematic diagram of a free path: O is the position of the (i-l) 
collision; M is the interface point of medium A and the clear region in the free path; M' is 
the interface point of medium B and the clear region in the free path; F is the position of 
the i collision. 






If b ext> x is not constant, one of the possible way to handle the free path is by a 
discrete method. We can divide the medium into boxes in which b ext x is a slow varying 
variable. Thus, b ext x can be seen as a constant in each box. Finally, / can be solved 

exactly the same way as equation (2.15). Unfortunately, calculating / discretely uses 
considerable computer time. In this research, we use the extinction coefficient at the (i- 
1) collision as an average extinction coefficient between the (/- 1) collision and the i 
collision. In another word, we use OF as the free path of OF. The drawback of this 
simplification is discussed in Chapter 4. 

As for the consideration of cloud shapes, it is actually an extreme case of b ext x 
variability. We define the cloudy region as b ext x * 0. and the clear region as b cxt j = 0. 

Thus, we can draw an imaginary interface between the cloudy and clear regions, which 
depicts the shape of clouds. Even though in the real world there is no such thing as an 
absolute interface to separate the clear and cloudy regions but a sharp decline of liquid 
water content (LWC) at the edge of clouds, we think this way of defining cloud shape 
serves as a satisfying approximation to study the effects of cloud geometry to the 
radiation field. 

Hence, when photons travel in the clear region, they would never collide with any 
cloud particles. Take the previous example of a photon traveling from medium A to 
medium B, and insert a clear region between medium A and medium B (Figure 2.3. b). 
The free path is just: 

/ = OF = OM + MM' + MF 
= OM (1 - > + _J_ i n (_L) +MM' 

k e xt.B k e xt,B 


. (2.16) 


Irregularity is the trade mark of cloud shapes; however, the general pattern of 
cloud shapes is of the main interest in our study. We can divide the cloud shapes present 
in the literature into three categories; (1) an isolated block (Davies, 1978; McKee and 
Cox, 1974); (2) organized and pattern-repeated cloud group (for instance, Blerkom, 
1971; and Welch and Wielicki, 1984); (3) stochastic cloud patterns (Zuev et al„ 1987). 
The stratocumulus clouds that we simulate in this study are well fitted by the category 2. 

In this study, our principles for handling different shapes of clouds are quite 
similar to Welch and Wielicki's research (1984); however, we add a few procedures as 
follows: 

1. Determine the smallest element of the alignment pattern called the unit cell 
(Figure 2.4). The highest point of the cloud is set to be Z^. 

2. Regardless of the variation of cloud height in this element, divide this element 
into boxes of equal size (Figure 2.4). 

3. Insert photons into each box at the center of each box at Z top . Calculate the 
location that photons would enter into the cloud as the starting point of the photon 
tracking process. 

4. Track photons. 

5. If a photon leaves the cloud, check if it would enter the adjacent cloud element 
or if it has already been reflected or transmitted. 

6. Calculate the free path of a photon that crosses the clear region. 

7. Apply the concept of periodic boundary conditions so that a photon that enters 
an adjacent unit cell is wrapped around to re-enter the original cloud unit cell (Figure 
2.4). 

8. Go back to step 4. 

9. After tracking all the photons, count the number of photons exiting from the 
top of each box (at Z top ). 


<7 



Figure (2.4) A schematic diagram for the unit cell, equal size boxes inside the unit 
cell, and how photons are wrapped to the unit cell when they go out to adjacent unit cell. 
In the diagram, the unit cell is uniformly divided into 16 photon collecting boxes. The V 
marks where incident photons enter the unit cell. A sample of photon wrapping is denoted 
by the thick line: O is the starting point; F is the real destination; F is the destination after 


wrapping. 



Figure 2.5 is a 


10. The reflectivity of each box is just (N out up / ^ ^ 

concise flow chart for the whole procedure. 

The bulk reflectivity of this particular cloud pattern is just the mean of the 
reflectivity of the equal-area boxes, or the total number of photons exiting from the top 
divided by the total number of photons inserted. Comparisons of our results with Welch 
and Wielicki's research (1984) are displayed in Figure (2.5) and (2.6). 

2.3.2 Statistics 

This leads us to the statistics question. What is the size of the photon collecting 
boxes and the number of incident photons per box needed to capture the 2-dimensional 
probability density function of reflecting photons, p(x,y; x',y’), which is itself a function 
of the location where photons are inserted? Here, (x,y) is the location where photons 
leave the top of the box, and (x',y') is the location where photons are inserted into the 
cloud element. 

To calculate the reflectivity of a homogeneous cloud deck, we just count the 
photons exiting from the cloud; we don't care where they come out. However, the 
location where photons exit is important in the case of inhomogeneous clouds. 

The analysis of the 2-D distribution of reflecting photons from homogeneous 
clouds provides a stepping stone for this problem. Figure (2.7.a), (2.7. b), and (2.7.c) 
respectively show the cumulative spectrums along the x axis, y axis, and the 2-D 
distribution of photons reflected when 2500 photons are incident on a homogeneous 
cloud deck at (0,0). Also, the width of the photon collecting boxes is 15 mean free path. 

Photons were inserted with a direction pointing toward negative x; hence, the 
cumulative spectrum along y axis is more symmetric than that along x axis. It is logical 
to think that if we insert enough photons, we can resolve the 2-D probability density 
function within satisfactory degree. Moreover, the smaller the collecting box size is, the 



The fraction cloud coverage 


Figure (2.5) Ratio of cubic field radiative fluxes to plane-parallel fluxes as a 

function of cloud cover: Welch and Wielicki (1984), this work. The incident 

zenith angle = 120°; aspect ratio of the cloud is 1/2, b ex , is 49 km* 1 .. 




The fraction cloud coverage 


.5 


Figure (2.6) Ratio of cloud field fluxes scaled to plane-parallel fluxes as a function 

of cloud cover for hemispheric clouds: Welch and Wielicki (1984); this 

work. The incident zenith angle = 120°; clouds are hemisphere shaped; b cxt is 73 km* 1 . 



Figure (2.7. a) The cumulative probability density function in the x direction for an 
infinite, homogeneous, non-absorbing cloud deck. The incident azimuth angle is 15°, the 
cloud optical depth is 16, g is 0.86. 



Figure (2.7.b) The cumulative probability density function in the y direction for an 
infinite, homogeneous, non-absorbing cloud deck. The incident azimuth angle is 15°, the 
cloud optical depth is 16, g is 0.86. 



Figure (2.7.c) The 2-D probability density function for an infinite, homogeneous, 
non-absorbing cloud deck. The incident azimuth angle is 15°, the cloud optical depth is 
16, g is 0.86. 


more incident photons are needed in order to get the detailed structure of the 2-D 
probability density distribution. 

For inhomogeneous clouds, every point at which photons are inserted into the 
cloud has its corresponding 2-D probability density function of reflecting photons, but 

the function is more skewed or asymmetric than that of the homogeneous cloud deck 
because of cloud geometry or varying b ext By overlapping the 2-D probability density 

O* P* 

functions, that is to say Ref p = J Jp(x,y;x',y')dxdy, we can obtain the reflectivity of 

each point. In addition, we want the size of the boxes large enough to capture the 
structure of the 2-D probability function and small enough to capture the geometry of the 
cloud. 

Now we can go back to explain the reason why we don't care where photons 
come out for the homogeneous cloud deck. Because the 2-D probability distribution is 
the same for every point in the infinite, homogeneous cloud, overlapping the probability 
distribution is equal to just counting the reflecting photons. From the above discussion, 
the Gaussian approximation to the binomial distribution is not applicable to describe the 
statistics of Ref b because Ref p is not a constant. Nor is it suitable for Ref p because 

Ref p = J J p(x,y; x’ , y' )dxdy . 


Chapter 3 


FIRE DATA 


The radiation and cloud microphysics data collected by the NCAR Electra flight 
2, which was carried out on June 30, 1987 during the FIRE marine stratocumulus 
observations, are of prime interest in this study. A GOES satellite image for flight 2 is 
shown in Figure (3.1), in which the flight path has been dotted. The cloud field can be 
categorized as a solid stratocumulus deck with occasional breaks. A lidar flight leg 
(referred to as leg 5 in the later text) with its corresponding turbulent flight leg (referred 
as leg 3 in the later text) are chosen for our simulation. 

Because of the differences in time and space of these two flight legs (Figure 3.2), 
we can not overlap the two data sets and obtain the correlation of the microphysics data 
and radiation data. However, some interesting information can be extracted from these 
data. Section 3.1 and 3.2 will describe the general cloud condition, and instruments used 
in leg 3 and leg 5. 

3.1. Lee 3 

Leg 3, in which the aircraft flew at a mean altitude of 790m above sea level, was 
carried out from GMT 20:10 (hounmin) to 20:20, which happens to be coincident with 
the time the satellite image (Figure 3.1) was taken. As it is shown from Figure 3.1 and 
from the summary of the cloud conditions in (Kloesel et al., 1988), the top of each 
individual cloud unit was solid, yet the cloud field as a whole was broken. 





Figure (3.1) The satellite image from GOES. 






The cloud microphysics data - the mean diameter (d), standard deviation (o), 
number density (N d ) and the liquid water content (LWC) - measured by the forward 
scattering spectrometer probe (FSSP) and by the 260X one-dimensional optical array 
probe are exploited. The FSSP, which has 15 channels, primarily sizes and counts water 
droplets from 0.5 pm to 45 pm in diameter, while the 260 sizes particles from 10 pm to 
620 pm in 10 pm increment. 

The definition of d and o follow basic statistics, d is simply the arithmetic 
average of all cloud droplet diameters: 

d=Xf,d, (3.1) 

i=l 


And c is calculated as: 


o = ^(Sf 1 d i 2 )-d 2 (3.2) 

And LWC is computed by: 


LWC = p „7Ef,d, ! 

6 i=i 


(3.3) 


, where fj = number density of particles in channel i, dj = diameter which channel i 

m 

represents, m = number of channels, and p w = the density of water. And N d = fj . 

i=l 


It should be noted that the LWC is not equal to 



Also, these definitions 


of d and a are different from the definitions of the geometric mean diameter and the 
geometric standard deviation of the log-normal distribution of cloud droplets, which is 
commonly used in calculating the extinction cross section, the asymmetry factor and the 



single scattering albedo in a standard Mie scattering code. See Section (4.1.1) for 
detailed illustrations. 

Figure (3.3.a) and (3.4) show d , N d and LWC measured by the FSSP and by the 
260X along leg 3. However, the x axis is intentionally reversed, because the flight 
directions of leg 3 and leg 5 were opposite (Figure 3.2). 

The N d and LWC in the Figure 3.3.a have the same tendency to increase from 
cloud edge to cloud center and decrease from cloud center to the opposite cloud edge, 
especially between 47 km and 39 km, and between 32 km and 27 km. Yet, d stays fairly 
constant inside of these two cloud blocks despite the existence of small fluctuations in its 
value. The other cloud blocks did not show apparent correlation among d , N d and LWC. 
The d and c are plotted in Figure (3.3.b). The a was magnified 10 times in order to 
show its correlation with d. Inside the two larger clouds (Figure 3.3.b), <7 does not 
change much. 

However, d , N d and LWC measured by the 260X (Figure 3.4), which is designed 
to measure the larger cloud droplets, did not demonstrate particular correlation in any of 
the cloud blocks. 

3.2 LegS 

Leg 5, in which the aircraft flew with a mean altitude of 1400m above the sea 
surface, or about 530m above the average cloud top height, was conducted from GMT 
21:09 to 21:20, approximately one hour after leg 3. The stratocumulus clouds of leg 5 
became more solid than the previous flight leg (Kloesel et. al., 1988). 

Figure 3.5 shows the equivalent brightness temperature, shortwave broad-band 
albedo and the cloud top height collected during this flight leg. The equivalent 
brightness temperature was transformed from the upwelling radiance measured by a 2°- 
field-of-view nadir radiometer within the atmospheric window region (10 - 12 fim). 



3 

8 

•S 

u 


•3* 

1 

o 

c 


§ 

s S 


CO 

2 

3 

W) 

E 


•s 

u 

s 

■s 

*0 


T3 

X 

O 

W5 


e ^ 

* i 

O U 

w *t3 

>» 
Xi 

1 
o 
c 
o 

T3 




£ 

o 


<D 

3 


0) 

5 

OX) 

c 

o 

T3 

cu 

00 

00 

uu 

0) 

X 

T3 

55 


I 

3 

j= 

3 

X) 


J2 

*D 

E 

3 

2 

u, 

3. 

X 

E 

3 

E 

O 

w 

* 

<L> 

•5 

bJQ 

3 


>> 

3 

a 

3 

X 

•o 

<D 

.£ 

73 

CO 

£ 

' — ✓ 

<3 

W) 

U* 

s 

o 

w 

o 



c 

o 

<N 

c3 

T3 

w 

X 

cn 


* * 

w> 


g S 

■§» g 

S) 3 

- 1 w 

| 5! 

3 £ 


o 


o 



along the NCAR Electra flight 2 leg 3 during FIRE. 




Brightness Temperature 



e 

w 

X 


1 

S3 

5/5 

o 

1 


c 

0) 

T3 

U 

•4-J 

<D 

T3 

E 

O 

*— » 
o 

o 

M 

XJ 

■o 

.5 

■S 

£ 

13 

*o 

o 

*■* 

*-» 


o 


-o 

c 


■s 

§ 

.S 

T3 

e 

•S 

i 

3 
«— > 
cd 

u 

5 

8 

; s. 

1 

3 

c3 

E 

t- 

8. 

WD 

S3 

60 

c 

E 

a 

C 

*C 

V5 

43 

60 

•c 

*o 

3 

T3 

V/~5 

s 

s 

43 


60 

60 

r 

jg 

•c 

43 


43 

60 

<N 

W3 


w 

O 

43 

43 

w 


60 

O 

o- 

o 

3 

c 

o 

■*-> 

Ctf 

T3 

1 *o 

is 


3 
1 3 

3 

1 

! ^ 

u 

TD 

i o 

a 

1 s 

3 

O 

W3 

i 

u 


/-N 



cr! 

'w' 

.s 

•+-T 

43 

2 

6 

C^J 

60 

*s 

3 

60 

o- 

43 

E 

<u 

43 

cu 

o 


w 

4 -» 


60 

T3 


C 

3 


o 

JO 


13 

13 


jtjgpH doi pnoio jo ‘opsqjv 


Shortwave broad-band albedo was calculated by dividing the upward shortwave flux 
measured by a downward-looking broad-band pyranometer by the downward shortwave 
flux measured by the similar instrument looking up. It has been verified that the 
downward flux was rather uniform along the flight path; therefore, there were no clouds 
necessary to be considered above the stratocumulus cloud deck. A downward looking 
lidar measured the cloud top height. 

As is shown in Figure 3.5, generally, the equivalent brightness temperature shows 
good positive correlation with the cloud top height. However, at the edge of the clouds, 
(holes in the cloud field are around 10-11 km and 37-38.5 km in the plot), the slope of 
the cloud top height is sharper than that of the equivalent brightness temperature. The 
explanation for this phenomenon is that lidar is a fairly sensitive detector of the existence 
of cloud droplets. The cloud top height is essentially the height at which the lidar detects 
a return signal from some volume containing cloud droplets. It is believed that there are 
still some cloud droplets scattering randomly in the breaks of clouds. Thus, using the 
cloud top height to infer the shape of cloud could underestimate the width of holes at the 
cloud top. 

The detected brightness temperature, however, relies much more on the cloud 
optical depth beneath the aircraft. For two clouds with the same LWC but with different 
number density, the cloud optical depth of the cloud with greater number density is 
larger. Yet, two clouds with the same number density but with different LWC, the one 
that contains more liquid water has larger optical depth. At the edge of clouds, the 
mixing of the moist air of clouds with the dry entrained air could result in a decrease of 
number density but an increase of the mean diameter. Both of these two changes will 
cause a decrease of the cloud optical depth at the edge of clouds. 

The most interesting feature of Figure 3.5 is that the albedo changes drastically 
from 0.1 to 0.65 while cloud top height and equivalent brightness temperature are fairly 
smooth except at the edge of clouds. The albedo maxima at 20 km and 30 km are 


coincident with the two small increments in cloud top height in the second cloud, 
whereas no obvious cloud top maximum is observed corresponding to the albedo 
maximum of the third cloud. It can also be recognized easily that the albedo declines 
almost linearly from the peak to the trough, or from the center of clouds to the edge of 
the clouds. 

The decrease of albedo at the cloud edge could result from the decrease of the 
optical depth at the cloud edge. Moreover, some photons which were originally inserted 
on clouds will leak out at vertical sides at the cloud edges causing additional decrease of 
albedo. The following chapter will discuss these two assumptions based on our 
simulation results. 



Chapter 4. 


Simulation Results 

Our motivation for the simulation originated from the plots of solar reflectivity 
and cloud top height along the flight path of FIRE Electra flight 2 Leg 5 (Figure 3.5). As 
it is shown on the figure, the cloud top height except at the edge of the cloud is rather 
uniform while the solar reflectivity shows a wave like pattern, which has its trough at the 
hole and its ridge at the cloud center. 

Two factors could contribute to this phenomenon: the geometry and the 
alignment pattern of clouds, and the variation of cloud optical properties inside the cloud. 
The goal of these simulations is to identify which of the two factors is more important in 
this case. 

4.1 Strategy 

Basically, two types of cloud field are used in our simulations of the visible 
radiation field. The cloud geometry and alignment pattern are considered in all of our 
Monte Carlo simulation experiments, but each experiment is with different degree of 
geometry effect. However, b ext remains constant in the first type of experiments, but it 
varies in space in a specified pattern in the second type. By comparing the result of each 
experiment with the observational data, we can obtain insight to the above cloud shape 
versus cloud optical properties problem. 

Meanwhile, the IR observational data can serve as an additional check to our 
simulations. We can calculate the IR radiation field corresponding to each simulation 
experiment and compare it with the observational data. Beside this, we can also check 



our assumptions of optical properties of the clouds with the cloud microphysical 
properties observed. 

In order to fulfill the comparisons mentioned above, a standard Mie scattering 
code, the IR hemispheric mean code, and a flux extrapolation scheme are introduced into 
this study. The following sub-sections briefly discuss these schemes. 


4.1.1 Mie s cattering code 


The volume extinction cross section o cxt (unit = [L 2 ]), the asymmetry factor g, 
and the single scattering albedo GJq of a specified log-normal distribution with geometric 
mean radius r g and geometric standard deviation a g are calculated by the Mie scattering 
code. The log-normal distribution of cloud droplets is assumed as follows: 


N(lnr)= vs5: expl i 




lna„ 


(4.1). 


The relation between r , the mean radius of cloud droplets (d/2 , d is defined in 
Section 3.1) and r g is r = r g exp[0.5(lna g ) 2 ], when the log-normal distribution of 

clouds is applicable. 

The theory of Mie scattering is not reviewed here. This sub-section only focuses 
on how we choose reasonable values for r g and o g to obtain satisfying values of a ext , g, 
andGJ 0 . 

Figure (4.1) series show plots of contours of c ext and g at 0.6 jim for different r g 
and c g . GJ 0 is not shown because it ranges from 0.999988 to 0.999999. The range of r g 
and G g in the plots covers a plausible range of r g and cr g in stratiform, non-precipitating 
water clouds. We can conclude that, for this type of cloud, the spectrum of cloud droplet 






GEOMETRIC STANDARD DEVIATION 


Figure (4. l.b) The contours of asymmetric factor at 0.6 ftmasa function of the 
geometric mean radius and the geometric standard deviation. Contour minimum: 0.852; 
maximum: 0.976; increment between two contours: 0.003. 



distribution does not affect the value of g and G3 0 much. The values GJ 0 = 1 and g = 0.86 
are fairly reasonable. 

However, the contours of a ext , g, and GJ 0 at 10 Jim are rather complicated. 
(Figure 4.2) We shall examine how we choose representative values of c ext , g, and GJ 0 
for our IR check in section 4.5. 

All the values of refractive index for water in this study are from Hale and Querry 

(1973). 

4.1.2 The Hemispheric M ean 2-Stream Technioue 

The Monte Carlo method is more complicated for IR radiation because the 
emission of photons also has to be considered. Therefore, the hemispheric mean 2- 
stream technique is selected to give us a first degree of approximation of the IR radiation 
field of a horizontally inhomogeneous cloud. The details of the hemispheric mean 
technique is described in Toon et al., 1989. A few approximations are made when the 
technique is applied to the calculation of IR field. 

First, we assumed that the temperature in the cloud deck is constant. This is 
actually a reasonable approximation. The value of a saturated adiabatic lapse rate at 
these temperatures is around 6°C in the lower troposphere; thus, for a cloud deck whose 
depth is only 200 to 300 meters, the change of temperature from cloud bottom to cloud 
top is only around 1.5 °C. 

Second, the ocean is assumed to be a blackbody. In another words, the 
absorptivity and emissivity of the surface are equal to 1 in our IR calculation. 

Then, we can write the equations as follows: 


F = kje* + rk 2 e' ET + C + 


(4.2), 









Figure (4.2.c) The contours of asymmetric factor at 10 Jim as a function of the 
geometric mean radius and the geometric standard deviation. Contour minimum: 0.888; 
maximum: 0.952; increment between two contours: 0.004. 



F = Fkje^ + k 2 e eT + C' 


(4.3), 


where F = upward flux, 

F = downward flux, 

C + = C- = 27rp 1 B(T c ), 

m = 1/2, 

T c = cloud temperature, 

B(T C ) = plank function at temperature T c , 
e = (7i 2 - 72 2 ) 1/2 » 
r = 72 / (7i + e), 

71 = 2 - 0J 0 (1 + g), 

7 2 = GJ 0 (1 - g). 

The boundary conditions are substituted into (4.1) and (4.2) to obtain kj and k 2 : 


F(t = 0) = 0, 


F (X = x c ) = 7tB(T s ), 

where x c is the IR optical depth of the cloud, and T s is the surface temperature. 

4.1.3 The Flux Extrapolation Scheme 

The purpose of this scheme is to extrapolate the reflectivity at the cloud top to 
some distance above the cloud deck. The observed albedo data was acquired when the 
aircraft was around 500 m above the cloud. We have to consider the extrapolation 
because the horizontal variability in reflected flux at the cloud top. 

In this scheme, first, we assumed that clouds are diffuse reflectors. That is to say 
the intensity reflected is isotropic. Therefore, F(x,y,z 0 ) = rt^x.y^). However, the 


intensity at z = z 0 + Az (where az is some incremental height above the cloud top) is not 
isotropic, I(x, y, z 0 + Az;0,<{>). The flux at z = + az is then equal to 

F(x,y,Zo + az) = JRx.y.Zp + az; ©,(|)) cos 0 sin 0 d0 d<J) (4.4). 

We have assumed that no scattering and absorption occur within the az atmosphere 

f 

column above the cloud. There are two ways to solve this integral. The first is more 
mathematical. 


Ffx.y.z, + az) = J J dx' dy (4.5). 


See Figure 4.3 for illustration. This is only a change of variables for the double integrals 

3(0, <b) . 30 d<{> 30 3<|) 

by using the Jacobian factor. The Jacobian — — is simply equal to — — . 

3(x' ,y) 3x' 3y 3y 3x 

Thus we can rewrite equation (4.5) as: 


Az 2 

F(x,y,Zo + az) = J J I(x',y,Zo) — — 2 , 2 dx' dy (4.6). 

((x -x) +y +z ) 

The second derivation contains more physics. We can rewrite (4.4) as: 

F(x, y^ + az) = J 1(0) cos © dO (4.7), 

because dO = sin 0 d0 d<|>. Since a horizontal surface element dA = dx' dy, on the z = Zq 

dA 

plane subtends a solid angle dO = — j 5 * 12 -, where dA ortho is the projection of dA 
orthogonal to the direction r, (4.7) becomes: 

F(x,y,Zo + Az) = J Itx'.y.Zo) cos © 

r 


(4.8). 


I (x,y,Zo+Az;0) 



Figure (4.3) A schematic diagram of the geometry of the radiance at cloud top 


and the radiance at Az above the cloud. 




Furthermore, because dA ortho = cos 0 dA, 


F(x,y, Zo + Az) = JJ«^^dx'dy 


(4.9). 


By substituting cos 0 = — , and r = ((x’-x) 2 + y 2 + Az 2 ) 1/2 , equation (4.9) eventually 

r 

becomes (4.6). 

Moreover, in our simulations, the intensity at the cloud top is constant with 
respect to y, because cloud geometry does not vary in y direction. (See section 4.2 for 
detailed descriptions.) Equation 4.6 can be further simplified into: 


F(x,y,z + Az) = J Kx'.Zq) {J 


Az 2 

[y 2 +((x'-x) 2 +Az) 2 ] 


dy) dx' 


I 2 * [(x ._ x) 2 +Az i]m 


dx’ 


— J R(x , ,Z 0 ) ^ T~rpr dx’ 

2 ^ [(x’-x) 2 +Az 2 ] 3/2 


= — E R(x i ,Z 0 ) ^ AX (4.10), 

2 ‘ ^ [(Xi-X) 2 +Az 2 ] 3/2 

where F 0 = incident solar flux on cloud, R(x i ,z 0 ) = the reflectivity at (x^Zq). Equation 
(4.10) is the equation we use to convert the reflectivity at the cloud top calculated by the 
Monte Carlo method into an upward flux at 500 m above the cloud top. 

The upper limit of Ax for the convergence of (4.10) has been carefully examined. 
We use an infinite homogeneous cloud deck of which the reflectivity is constant in space 


as our first-step approach to the convergence problem. It is found that when az = 500m, 
200m is the upper limit of ax. Therefore, the size of the photon collecting boxes could 
not exceed 200m in width in our simulations. 

Az 2 

From another perspective, = 7-77 ax serves as a weighting factor of 

2[(Xj -x) +Az ] ' 

the reflectivity at R(Xj,Zo) for the calculation of F(x,y,z„ + Az). Let's see the situation 
when Az = 500 m and ax = 200 m. When (x r x) is equal to zero, the weighting factor is 
0.2; when (x r x) is equal to 1 km, the weighting factor drops to 0.018; and when (x r x) is 
equal to 5 km, the weighting factor further decreases into 0.0002. Furthermore, we have 
estimated that if the intensity were isotropic, about 90 % of the measured flux would be 
contributed by reflected intensities from cloud top within a width of 1.1 km around the 
aircraft which is at a height of 500 m above cloud top. 


CT t \ fill? f5i iTiT* 


From the satellite image (Figure 3.1), we can see that the clouds along the path 
that the airplane flew are of the cloud-street type. Thus, in our simulations, the x axis is 
defined as the flight path. We assumed that the clouds extend to infinity in the y 
direction, and the geometry and cloud optical properties are only functions of x. That is 
to say, the unit cell (defined in Section 2.3.1) is an infinite cloud strip extending along y 
axis. As a result, the reflectivity at the cloud top remains invariant in the y direction. 
Because of this simplification, we need only collect outgoing photons along the x axis 
when they leave the cloud. 

Table 4. 1 lists all the constants necessary to be specified in order to execute the 
simulations and the values that we choose. 

The incident solar zenith angle 0 O and the incident solar azimuth <{> 0 are calculated 
by the following equations (Spencer, 1971): 





Descriptions 

Values 

Incident angles 

€>o 

Incident zenithal angle 

165° 


4>o.. 

Incident azimuthal angle 

90° or 0° 

Cloud geometry 

H 

Cloud depth 

260 m 


w h 

Width of the hole at cloud bottom 

1km 


W c 

Width of the cloud at cloud bottom 

26.1km 


1 

Asymmetry factor 

0.86 



Single scattering albedo 

1 


^ext, X^0.6 \im 

Volume extinction coefficient 

0.06154 m- 1 * 


h 

Mean free path 

16.25 m 

Box width 

Ax 

Width of a photon collecting box 

200 m 

Surface albedo 

A, 


0 


Table (4.1) The list of all the constants in the simulations and their values . *In 
experiment 3 and 4, the volume extinction coefficient decreases linearly, or quadratically 
from 0.06154 at the cloud center to 0 at cloud edges. 




















cos 0 O = sin 0 sin 5 + cos 0 cos 8 cos h 


(4.11), and 


sin <(> 0 = -cos 8 sin h / sin O 0 (4. 12). 

In the above two formulae, 0 is latitude. 

8 (in radian) = 0.006918 + 0.399923 cos <p 0 + 0.070257 sin <p 0 - 0.006758 cos 2cp 0 + 
0.000907 sin <p 0 - 0.002697 cos <p 0 + 0.001480 sin 3<p 0 (4. 1 3) 

(<p 0 =2ixd n / 365, d n is the day number ranging form 0 on January 1 to 364 on December 
31), and h is the hour angle, which is symmetrical about solar noon in terms of the local 
apparent time. 

local apparent time = clock time + longitude correction + 0.000075 + 0.001868 cos(<p 0 ) - 
0.032077 sin (cp 0 ) - 0.014615 cos (2cp 0 ) - 0.04849 sin (2<p 0 ) (4.14) 

However, we did not use the <j> = 295° (0° is south, and the angle rotates 
counterclockwise), calculated according to Spencer’s formulae. The reason is that there 
is too much uncertainty embedded in calculating the incident azimuth angle given our 
arbitrary coordinate system. Instead of using an uncertain azimuth angle, we decide to 
input photons at <|> 0 = 90° and 0° (0° points to the positive x direction, 90° points the 
positive y direction and so on) for the incident azimuth angle in our simulations. In the 
case of 90° azimuth angle, the vertical cloud side effect should be the smallest, and at 0° 
it should be the largest. 

H, W h and W c are obtained by averaging leg 5 data. However, the ceiling height 
was determined by another earlier flight. 

Section 4.1.1 has already discussed the values for g and GJ 0 . Instead of averaging 
the microphysics data to acquire b ext ^ at visible wavelengths, the b ext ^_o 6)am was 

estimated from the reflectivity data. From the Monte Carlo code for homogeneous 



clouds, we found that the optical depth is around 16 for clouds of reflectivity 0.55, which 
is the mean value of reflectivity around cloud center (Figure 3.5). We then divided cloud 
depth (260 m) by the optical depth 16 to get the mean free path. Discussions about this 
simplification is in Section 4.3. 

The upper limit for the convergence of our flux extrapolation scheme, 200m, is 
used as the width of the photon collecting boxes. The ocean is assumed to be black in 
the visible wavelengths. In fact, this assumption is reasonable (Figure 4.4) since the 
zenith angle is small in this case. 

Table (4.2) summarizes the four simulations we conduct. Experiment 1, in which 
the geometry effect is minimum and b ext is constant, serves as the control mode of the 
simulations. In experiment 2, the cloud geometry effect is a maximum. In this run, the 
cloud geometry is most similar to the cloud top height data. Experiment 3 and 4 
investigate the effects of a variation of b ext on reflectivity with small cloud geometry 

effect. Figure (4.5) displays the three kinds of cloud shapes used here in comparison 

* 

with the cloud top data of leg 5. 

In experiment 3 and 4, because of our simplification of tracking photons through 
a b„ t varying medium (Section 2.3.1), one potential problem rises when photons collide 
with cloud droplets near cloud edges. At the edge of cloud, the mean free path (1 / b Mt ) 
will increase drastically; therefore, we will overestimate free paths (Equation 2.12). 

4.3 Results 

Figure (4.6) shows the results of the four simulations. The azimuth direction of 
the incident photons points to the positive x direction in each case. 

As it is shown on the plot, the results of simulation 1 and simulation 2 are similar, 
even though the reflectivity of simulation 1 varies slightly sharper in the vicinity of the 
cloud edge. This is not a surprising result, because, from Figure (2.7), the cumulative 
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Figure (4.4) The albedo of the ocean surface as a function of the solar zenith angle 

and the vertical transmissivity of the atmosphere: , clear (transmissivity 0.6-0.65); 

, thin overcast (transmissivity 0.35 -0.4), , computed for a parallel beam and a 


smooth surface. Courtesy (Kondratyev, 1969) 



# 


Descriptions 

kr.Tf 

Descriptions 

1 

Rectangular 

Rectangular cloud 
blocks 

Constant 


2 

Trapezoidal 

Trapozoidal cloud 
blocks with a width 
of 26. 1 km at bottom 
and a width of 22.1 
km at top 

Constant 


3 

With 

rounded edge 

Cloud edge is 
rounded with a radius 
of 260 m 

Linear 

variation 

b ext decreases linearly 
from 0.06154 m 1 to 0 
from cloud center to 
cloud edge 

4 

With 

rounded edge 

The same as 
experiment 3 

Quadratic 

variation 

b ext decreases 
quadratically from 
0.06154 nr 1 to 0 from 
cloud center to cloud 
edge 


Table (4.2) Descriptions of 4 experiments. 


















dashed line for experiment 2; the dotted 



probability density diagram in the x direction, in which the optical depth is 16, and the 
incident azimuth angle is 15°, shows that the probability of photons reflecting from the 
cloud at a distance 75 times the mean free path is very small. Thus, the geometry effect 
could only possibly affect the radiation field about 1 km from the starting and the ending 
points of the geometry variation. For example, for experiment 1, the geometry only 
affects the reflectivity between 36.1-39.1 km, for experiment 2, between 34.1-41.1 km. 

In contrast, the results of simulation 3 and 4 show more gradual change of 
reflectivity from the cloud center to edge. Except near the center of the clouds where the 
reflectivity of experiments 3 and 4 are almost the same, the reflectivity of experiment 3 is 
less than that of experiment 4. This is consistent with the optical depth variation of the 
two experiments. (Figure 4.7) 

Similar experiments were conducted to investigate the azimuth angle effects. 
This time, the incident azimuth angle was set to be 90°; that is to say, the incident 
azimuth angle is directed along the positive y direction. Figure (4.8) is a series of 
comparisons of two extreme incident azimuth angles for each of the four experiments. 
From figure 4.8, there are no obvious differences between the 90° and 180° cases. The 
original conjecture of these experiments was that the reflectivity of 180° incident azimuth 
angle would be larger than that of 90° cases at the vicinity of cloud, because the vertical 
sides of clouds would reflect more photons for 180° cases. However, we can not see 
obvious differences at the vicinity of clouds. Possible reasons are the incident zenith 
angle is too close to that of solar noon, and the fractional cloud coverage is near unity. 

Since the difference between 90° and 0° cases is not obvious, only the 0° cases 
are used in the discussions below. We now apply the flux extrapolation scheme to the 
raw results. Figure (4.9) is a series of comparisons between the cloud-top results and the 
extrapolated results. From these figures, we can conclude that the flux extrapolation 
scheme smoothes the variability of the Monte Carlo method. Also, the reflectivity of the 
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Figure (4.7) A plot of optical depth for four experiments. The dashed line denotes 
experiment 1; the dotted line for experiment 2; the solid line for experiment 3; and the 
dashed-dotted line for experiment 4. 
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Figure (4.9.c) Albedo at the cloud top before extrapolation (the light line) and 
albedo after extrapolation (the heavy line). Results from experiment 3. 
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hole increases after extrapolation due to side effects. This is consistent with our 
discussion of the flux extrapolation scheme in Section (4.1.3). 

We can plot these modified results with our observational data (Figure 4.10). The 
observed albedo matches fairly well with the experiment 4 in the left part of the cloud. 
At the hole, both experiment 3 and 4 underestimate the reflectivity but only with an error 
about 0.05, which is close the 5 % sea surface reflectivity for small azimuth angles in 
Figure (4.4). At the right cloud block, the observed albedo first agrees with the 
reflectivity calculated by experiment 4 and then it becomes smaller than that of 
experiment 3 with a discrepancy about 0.05 in maximum. Both experiment 1 and 
experiment 2 are unable to create a gradual change of reflectivity from cloud center to 
cloud edge. They only match the observed data well at the hole. Thus, this implies the 
geometry effect alone is not sufficient to result in the albedo change of observation. 

Now, the IR data are used to check our results. First of all, we have to determine 
the g, GJ 0 c ex t for IR at 10 |im. The log-normal distribution of cloud droplets 
assumed in the Mie code is checked from data (Figure 4.11). We found that the log- 
normal distribution assumption is good when the cloud is solid, or when the observed 
data is at the center of clouds. The cloud microphysical data is then used to obtain the 
geometric mean radius and geometric standard deviation for those sets satisfying the log- 
normal distribution. The points with "+" symbols in Figure (4.12) are sets of the r g and 

a g calculated from observational data. 

In section (4.3), we have mentioned that the value of the extinction coefficient in 
the visible is set to be 0.0654 nr 1 . However, from the observational data, the number 
density of cloud droplets is around 30 to 55 per cm 3 . Even if the maximum 55 per cm 3 is 
input for the calculation of extinction cross section: 


extinction coefficient = volume extinction cross section * number density (4. 1 1), 





Figure (4.11) The scatter plot of the standard deviation calculated from measured 
LWC and d by FSSP using the log-normal cloud droplet distribution assumption (y axis) 
and the standard deviation measured by FSSP (x axis). The straight line denotes cases 
that perfect matches of the two standard deviations take place. 






Figure (4.12) The scatter plot of the observed geometric mean radius and 
geometric standard deviation (denoted by V) which satisfy the log-normal cloud droplets 
distribution during leg 3. The units for both the geometric mean radius and geometric 
standard deviation are in pm. Meanwhile, the contours of the extinction cross section (p 
m 2 ) at 0.6 pm are also shown. 




the value of extinction cross section is about l.OxHT 9 m 2 , which does not fit our cloud 
droplet distribution data (Figure. 4.12). 

Yet, it is found that the CCN (cloud condensation nuclei) concentration in the 
boundary layer was around 70 to 150 per cm 3 (Albrecht, personal communication). If 
we choose 75 per cm 3 and substitute into equation (4.11), the extinction cross section 
matches well with the cloud droplet spectrum observed. The cause of the discrepancy 
between the observed cloud droplet number density and the CCN number concentration 

may be that the FSSP does not measure all droplets. 

From this procedure, r g = 9.5 pm and a g =1.4 pm are used to obtain the IR 

extinction cross section, IR asymmetry factor, and IR single scattering albedo. The 
values of these radiative properties are 900 pm 2 , 0.93, and 0.68 respectively. 

The 2-stream hemispheric mean scheme is then used to calculate the brightness 
temperature of these four types of clouds. The results are shown in Figure (4.13). 
Again, the geometry alone could not match the measured brightness temperature. 
However, when the extinction coefficient decrease linearly or quadratically from the 
cloud center to the cloud edge, the calculated brightness temperature are close to what is 
observed. 

In our experiments 3 and 4, the extinction coefficient varies. However, the 
asymmetry factor and single scattering albedo were held constant. Even though that the 
asymmetry factor and single scattering albedo remain rather constant for common cloud 
droplet spectrum for visible radiation, the variation of the asymmetry factor and single 
scattering albedo in the IR is not negligible. By this, in our simulation, the cloud droplet 
distribution is assumed to be invariant in clouds, while the number density changes from 
cloud center to cloud edge. In fact, if we go back to Figure (3.3.2), it is clear that 
between 47 km to 39 km and between 32 km to 27 km of our arbitrary coordinate that 
the r g and G g stay rather constant while number density increases from cloud edge to 
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x (m) 

Figure (4.13) The observed brightness temperature and the brightness temperature 

for the four experiments (°C) along the flight path (20 km to 60 km). The observed 
brightness temperature has distinct small fluctuations; the other four lines are experiment 
1, 2, 4 and 3 from the center of the hole inside out 



cloud center and decreases again near the cloud edge. However, because the time and 
space difference of leg 3 and leg 5, direct comparison is impossible. 

One of the potential problem in these simple simulations is that the vertical 
variation of cloud microphysical properties is not considered. It is suggested that the 
vertical variation of cloud optical properties could be included and the effective particle 
radius of clouds (Nakajima and King, 1990) could also be applied in the further 
simulation of the radiation field of broken cloud field. 

It should also be emphasized that, in this study, it is not the detailed quantities in 
which we are interested. But after analyzing of our results, we can conclude that for the 
case studied in these simulations, cloud geometry alone could not explain the drastic 
variation of reflectivity detected from the airplane, nor the brightness temperature 
variations. 

Further studies of broken clouds would be valuable especially in collecting cloud 
microphysical data; for instance, mean radius, standard deviation and number density as 
a function of relative position in clouds. Therefore, a better understanding of 
stratocumulus clouds could be achieved, because these simulations infer the drastic 
differences of radiation field between plane-parallel clouds and its corresponding broken 
clouds in which the cloud optical properties change from cloud center to cloud edge. 



Chapter 5 


CONCLUSION 


The Monte Carlo method was applied to simulate the visible radiation field over a 
broken cloud field. Four experiments have been conducted to reproduce the reflectivity 
and the brightness temperature measured in FIRE Electra flight 2 leg 5. One most 
significant feature in this flight is that the albedo varies almost linearly from 0.6 at cloud 
center to 0.05 at the holes between clouds, while the brightness temperature and the 
cloud top height remain unchanged through out the clouds. 

The results suggest that the cloud geometry alone could not count for the change 
of albedo and the brightness temperature from cloud center to cloud edge. Both 
experiment 1 and 2, in which only the cloud shape effect is considered, fail to mimic the 
observed data. In contrast, both experiment 3 and 4, in which both the cloud effect and 
the variation of optical properties are considered, get consistent results with the measured 
data. It is the change of cloud optical properties inside the clouds more possible to be the 
main cause of what has been observed. 

Assumptions and approximations in the simulations require further careful 
studies, especially for the building of a more accurate and fast way to generate free paths 
in a inhomogeneous medium. 

In experiment 3 and 4, we vary the extinction coefficient linearly or quadratically 
from the cloud center the cloud edges. Because the extinction coefficient is a product of 
the cloud droplet number density and the volume extinction cross section, extra caution is 
needed to infer our results to the variation of microphysical properties inside the clouds. 



In Section (3.2) we have demonstrated that the observed mean radius and its 
standard deviation in Flight 2 Leg 3 are rather constant in comparison to cloud droplet 
number density and LWC. Also, results of experiment 3 and 4, in which we assumed 
that the cloud droplet spectrum remained unchanged through out the cloud but only the 
cloud number density changed, are similar with the albedo and brightness temperature 
observed. 

Conventionally, we think the cell structure of stratocumulus clouds is caused by 
the rising motion at cloud center and sinking motion at cloud edges. The cloud droplets 
will grow in the ascending air; however, when they reach the cloud top, they will mix 
with the entrained dry and warm air which are from the inversion above the cloud and 
start to evaporate and sink. The adiabatic warming in the descending motion will further 
contribute to the evaporation of cloud droplets. Smaller cloud droplets will soon vanish 
in the descending branch of vertical motion; meanwhile, a shift of cloud droplet spectrum 
toward larger cloud droplets is expected. 

Yet, how does this argument link with the observed microphysical data and the 
simulation results? Why is the change of cloud droplet spectrum not observed? Why do 
we get good results when the change of cloud spectrum is not considered? Does that 
mean the shift of cloud droplet spectrum is small compared to the evaporation of smaller 
cloud droplets? More studies about the microphysical structure of stratocumulus clouds 
are needed to answer these questions. 

It is also conjectured that whether the case studied here represents a typical case 
for stratocumulus clouds or not. First, the case we studied was at around solar noon, the 
cloud vertical side effect was actually at the minimum situation. Second, the horizontal 
scale of the breaks of the clouds is much smaller than that of clouds. Third, the clouds 
are approximately semi-infinite in geometry, because their horizontal scale (~25 km) 
overwhelms the effective width (~1 km) of the 2-D probability distribution of its 
homogeneous counterpart. Indeed, for those stratocumulus clouds that the fractional 



coverage of cloudy area is smaller and that the size of the cloud is comparable to the 
effective width of their 2-D probability distribution of reflected photons, the cloud 
geometry effect could be much more important. However, if this case is representative, 
not only the effect of cloud geometry but also the variation of cloud optical properties 
must be considered in the climate models. 
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